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ABSTRACT 

In Newtonian gravity the final states of cold dissipationless collapses are characterized 
by several structural and dynamical properties remarkably similar to those of observed 
elliptical galaxies. Are these properties a peculiarity of the Newtonian force or a more 
general feature of long-range forces? We study this problem by means of N— body 
simulations of dissipationless collapse of systems of particles interacting via additive 
r~ a forces. We find that most of the results holding in Newtonian gravity are also 
valid for a ^ 2. In particular the end products are triaxial and never flatter than an E7 
system, their surface density profiles are well described by the Sersic law, the global 
density slope-anisotropy inequality is obeyed, the differential energy distribution is 
an exponential over a large range of energies (for a ^ 1), and the pseudo phase-space 
density is a power law of radius. In addition, we show that the process of virialization 
takes longer (in units of the system's dynamical time) for decreasing values of a, and 
becomes infinite for a = — 1 (the harmonic oscillator). This is in agreement with the 
results of deep-MOND collapses (qualitatively corresponding to a = 1) and it is due to 
the fact the force becomes more and more similar to the a = — 1 case, where as well 
known no relaxation can happen and the system oscillates forever. 

Key words: gravitation - stellar dynamics - galaxies: kinematics and dynamics — 
methods: numerical 



1 INTRODUCTION 

One of the most striking properties of elliptical galaxies is 
the remarkable quasi-homology of their surface brightness 



~J „„„ ~~ ~ 

alization of the de Vaucouleurs 7? 1 / 4 model (see e.g. 


Caon 


et al.| 1993| |Andredakis et al. 1995| Courteau et al. 


1996, 


Graham & Colless 1997, Prugniel & Simien|1997[ Graham 


1998 


|Trujillo et al.|2001, Bcrtin et al.|2002| see also|Ciotti 


2009 


and references therein). Albeit minor (but important) 



departures from the Sersic model are common, overall the 
profiles on large scale are very well represented by the 
Sersic law. What is the origin of such regularity? A-body 
numerical simulations revealed that cold dissipationless and 
collisionless collapses lead to virialized e nd-states described 
almost perfectly by R 1 ^ profiles (e.g. 



van Albada 



1982 



Londrillo et al. 19911. More recently it has been shown 
that collapses in pre-existing dark matter halos are also 
well described by the Sersic profile, with a wide range of 
values of the Sersic index (Nipoti et al. 2006ab, hereafter 
N06ab). In addition, it is also known that the Sersic 
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family is characterized by an exponential differential energy 
distribution, over a large range of accessible energies (e.g. 
Binney|1982] |Ciotti|199l"j ). These results can be understood 
in terms of the physics of violent relaxation in collisionless 
collapses (e.g. |Lynden-Bell||1967[ |Bertin fc Stiavelli||1984[ 



Bertin fc Trenti| |2003[ |Trenti fc Bertin| |2005[ |Trenti et] 
al. 20051. Finally, it has been proved analytically that in 



Newtonian gravity a large class of spherically symmetric 
equilibrium systems are characterized by the so-called Global 
Density Slope- Anisotropy Inequality (hereafter GDSAI, see 
Ciott i & Morg anti 2010ab; |van Hese et aLl|2011[ |An et al.| 



2012 see also An & Evans 20061, a constraint between 



their anisotropy and density profiles. Numerical simulations 
suggest that the GDSAI may be a much more general 
result, holding true also for the final states of dissipationless 
collapses (see e.g. Hansen fc Moore]|2006 1. 



Due to the relevance of these results for the under- 
standing of the process of collisionless relaxation, a natural 
question arises about their apparent universality. In particu- 
lar, are the Sersic law, the associated differential exponential 
energy distribution and the GDSAI peculiar features of New- 
tonian gravity or are they more general properties of the 
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virialized final states of TV-body collapses in which the par- 
ticles interact with long-range forces? Preliminary results 
seem to support the second possibility. For example, it is 
known that the end-products of cold collapses in Modified 



theorem on spherical harmonics, so that in principle a 



Newtonian Dynamics (MOND, |Bekenstein & Milgrom|l984 \ 
also produce final systems described remarkably well by the 
Sersic law ( Nipoti et al.|2 007a hereafter N07a; and Ciotti et 
aL]|2007[ ). However, the TV-body MOND simulations have 
also shown that the oscillations leading to relaxation last 
more (in units of the dynamical time of the system) than in 
the equivalent Newtonian system (see e.g. N07a; Nipoti et 
aT|2007b| ). Moreover [Barber et aL] ( |2012| ) found indications 
that the GDSAI may be a common property of MONDian 
virialized systems. We recall that the force law in the MOND 
weak field limit (or deep-MOND, hereafter dMOND) regime 
is qualitatively similar to a force decreasing with distance as 
1/r. 

Additional indications in this direction come from the 
preliminary analysis of Di Cintio (2009, hereafter DC09) 
and Di Cintio & Ciotti (2011, hereafter DCC11), who 
investigated the relaxation of a system of spherical shells 
interacting via a long-range r~ a force law, in analogy with 



similar studies performed in Newtonian gravity ( Henon 



Takizawa & Inagaki 1997 Youngkins & Miller 2000 1 and 



1964 



MOND (Sanders 1998, 2008; Malekjani et al. 2009, 2012). 
In DCC11 we focused on the time evolution of the virial 
ratio and of the differential energy distribution. Among the 
main results, we confirmed the expectation that the process 
of relaxation, independently of the value of a (with the 
exception of a = —1), consists in a first phase of violent 
collapse, followed by a longer, gentle phase of dynamical 
mixing. Remarkably, small values of a correspond to larger 
and long-lasting virial oscillations, confirming the dMOND 
results (N07a). However, the final states of shell systems 
are only poorly described by an exponential differential 
energy distribution. This is not surprising since the 
enforced spherical symmetry reduces the number of degrees 
of freedom available for energy exchanges during virialization. 

Prompted by these preliminary results, here we explore 
further the problem following the collapse and virialization 
of fully three-dimensional TV-body systems of particles 
interacting with radial forces proportional to a power-law 
of the mutual separation, r~ a . This approach is not new, 
and here we recall the study of quasistationary states 
( |Gabrielli et al.|2010[ [Marcos et al.|2012| and of Yukawa-like 
gravity (Moffat & Sokolov 1996 and references therein; 



see also Brandao & de Araujo 20121. For the simulations 



we developed a direct TV-body code, exploiting the force 
additivity. Note that for general forces, more sophisticated 
methods, based on the expansion in orthogonal functions 
of the potential, are not available, as the analogue of 
the Poisson equation does not exist. In our case, the 
considered forces even though additive, are described by 
non-local operators, i.e. the density at a given point can 
not be expressed as a simple differential operator of the 
potential at that point. In fact, the potentials associated 
with r~ a forces are the well known Riesz potentials, and 
the analogue of the Poisson equation involves the so called 
fractional Laplacian (e.g. Stein|1970 1. Remarkably, for this 
specific cases, the potential can be expressed in terms of 
Gegenbauer polynomials in turn expressible with an addition 



multipole based Treecode can be realized (see Srinirasan 
et al.|2005| ). Note that MOND is a non-linear but local theory. 



The paper is organized as follows. In Section 2 we in- 
troduce the most important integral identities that will be 
used to study the results of the simulations, while in Section 
3 the numerical code and the set-up of the initial conditions 
are presented. In Section 4 the virialization process and the 
structure and the dynamical properties of the virialized final 
states are presented and discussed as a function of a. The 
main results are finally summarized in Section 5. 



2 SETTING THE STAGE 

We integrate numerically the equations of motion for an 
initially spherical system consisting of TV particles of iden- 
tical mass m, mutually interacting with central long-range 
forces, obeying the superposition principle. In particular the 
acceleration at r; due to a particle of mass nrij at Vj is 



iji = —Grnj 



r; - r i 



(1) 



where | | is the standard Euclidean norm and G, is the force 
constant. The associated potential is 



ji = Gmj 



hi 



1 • 



a / 1, 

!=1. 



(2) 



with &ji = 



-X7i(j>ji. Note that for a = 2 we recover the 
Newtonian gravity, for a = — 1 Hooke's harmonic force, 
and finally a = 1 is the additiv^] analogue of the dMOND 
regime in which the MOND force behaves qualitatively as 
1/r. Some authors (e.g. Chavanis 2008 , [Bouchet et al.|2010 1 
introduce the nomenclature of weak long-range interactions 
and strong long-range interactions depending on whether 
the force vanishes or diverges for r — > +co; here we label 
them as gravity-like (a > 0) or harmonic oscillator-like 
(a ^ 0). The forces in eq.(l) can be also divided in two 
families depending on the confining nature of their potential, 
separated by the case a = 1, when the potential diverges 
for zero and infinite separation. For a 1 the total energy 
of a particle may be positive or negative, while it is always 
positive for a < 1 (having assumed zero potential energy 
for a system collapsed at the origin). Particles can escape 
only from systems with a > 1, while bound particles have 
negative energies. 

In order to keep track of the process of virialization, for 
each simulation we follow the evolution of the virial ratio 



1 



2K 

W\ 

where 

N 



rriiVi 



(3) 



(4) 



1 Recently, |Milgrom| j201u) proposed a quasi-linear formulation of 
MOND called QuMOND. We stress that the a = 1 case studied 
here is not QuMOND. 
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Figure 1. Newtonian (a = 2) tests for Plummer initial conditions with r/o 0.5. Left: projected density profiles for the final states 
at 50 t* (dots), and their best-fit Sersic profiles (lines); residuals are also shown. Right: normalized differential energy distribution 
n* (-E) — 7i(£/) j N for the final states. The black solid line represents n*(E') for the initial condition with 7?o = 0; note that E m [ n depends 
on the specific realization. 



is the total kinetic energy of the system, and 



|H]) one obtain^] 



W = - I p(x) < x, V0 > d 3 x = ^2 nii < Xi,a.ji > (5) 



is the virial function, so that the virial theorem reads 2K = 
—W. Note that the condition j 7^ i in eq. |5| is required for 
q > 0, however this condition can be extended without loss 
of generality also to a ^ 0, due to the vanishing of the self 
force. Note also that convergence of the integral in eq.|5| 
requires a < 3. We recall that in general W is not the total 
potential energy 



U 



If 1 

2 / p(x)(/>(x)d 3 x = - 2^ rrncpji, 



(6) 



but for a 7^ 1 it is proportional to it, being 



W = {<x- l)U. 



(7) 



Note that U — when the system is dispersed at infinity 
for a > 1, while U = when the system is collapsed at the 
origin for a < 1. In the a — 1 case, the reference state must 
be fixed with the particles at finite (but non zero) separation. 
The case a = 1, corresponding to the logarithmic potential 
(i.e. to dMOND-like force), is peculiar, as from eqs. |TJ and 



W 



G 



rriimj, 



(8) 



i.e. W remains constant during the virialization. Remarkably, 
it can be shown analytically that the virial function is time- 
independent also in dMOND, even though the field equation 
is not linear and the force is in general neither radial nor 
strictly proportional to 1/r, (N07a; for the special case of 
spherical systems see Gerhard fc Spergel|1992 \. As the total 
energy E = K + U is conserved for all values of a, oscillations 
in K are always associated with oscillations in U, so that for 
a 1 eq. |7| and the Lagrange-Jacobi identity I — 2(2^+1^) 
(e.g. Ciotti| 2000 1 show that the time dependence of the 



moment of inertia of the system is due to the combined 
effects of K and U. For a = 1 instead only the kinetic energy 
changes during the virialization. Finally, for a dissipationless 
collapse, starting from cold initial conditions (K^ = 0), it is 
easy to prove that the value of the equilibrium wirial function 
(Wfin) is related to the initial potential energy Ui n by 



W fin = 



2(a ■ 



3-c 



(9) 



where A"| n is the asymptotic energy of the possible escapers. 
In the identity above it is assumed that the escapers are fully 
dispersed, i.e., that their gravitational energy is small, and 

2 Note that in the case of a continuous density distribution with 
o = l and total mass M, W = —GM 2 /2, and this is not the limit 
of eq.(8) for N -s> 00 and m = M/N. 
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Figure 2. Final axial ratios as function of the force exponent a, 
for Hernquist and Plummer initial conditions with different values 
of t]q. Note how for each value of 770 the results are almost identi- 
cal independently of the initial density profile. In the harmonic 
oscillator case (a = — 1) the systems retain their spherical shape. 



that Wan is the virial function of the remnant. Not also that 
in the cases with a < 1 no escape is possible, so that eq.Q 
holds rigorously with K^ n = (provided equilibrium could 
be attained), and Wfl n refers to the whole system. We used 
eq.|9| as a test for the simulations. 



of the system (see also DCC11); the natural velocity scale 
becomes 



r. 



(11) 



Setting x = r/r* and r = t/t„, the dimensionless equations 
of motion for the particle i become 



"x, 



dr 2 



2 

"N 



E 



X, 



(12) 



We note that, as all the presented results are scaled to the 
dynamical time t* , the specific value of G does not affect the 
conclusions. The code used for the simulations is a direct code, 
with the consequent limitations on the number of particles 
that can be used: as a rule we adopt TV = 25000. In the case 
of gravity- like forces (a > 0), the divergence of the force 
when the interparticle separation tends to zero is cured with 
the introduction of the softening length e (e.g. Dehnen|2001 1, 
so that eqs.|l|-([2| are replaced by their softened expressions 



soft /~* 

a,, = — Gm, 



Aj° = Gnij x 



r, 



Ti - r, 



+ e 2 )" 



1 - a 
In a/ 1 1 r . ; - r j 1 1 : 



+ e 2 



a = 1. 



(13) 



(14) 



The optimal e is chosen as follows: the density profile of the 
initial conditions is divided in spherical shells of radius Ri 
and thickness SRi, and the minimum inter-particle distance 
di within each shell is computed. A value e, is obtained 
by comparing the acceleration a Boit (di,€i) of a pair with 
separation di with af° ft (e 4 ), the acceleration of a random 
particle of the shell due to all the other particles of the system. 
The value of is chosen so that a'° ft (ei) = a sott (di, £»). As 
optimal e we take the maximum e^. We verified that with 
such choice, independently of a, the softened acceleration 
of a particle at large distance from the centre of mass of 
the system differs by less than 0.01% from the non-softened 
acceleration. The equations of motion are integrated using a 
standard second order leapfrog method both with constant 
and adaptive timestep At. In the second case At = min(Ati), 
where AU — mm(Atu, At?i, At$i), and 



At u 



|Ari| 
lAvil 



At 2l 



x»|jAri| 
IIAJdl 



AU 



mi] | Ari 
IAEA 



Ar,, Avi, AEi and AJi are the variation of position, velocity, 
energy and angular momentum of the ith particle in the 
previous timesteps. 



3 THE SIMULATIONS 



3.2 Initial conditions 



3.1 The N -body code 

In order to compare the process of virialization of systems 
with different values of a, we introduce the time-scale i* 
from the relation 

GMtj 
2r? +1 



1, 



(10) 



where r, (the length scale) is the half-mass radius of the 
density distribution at t = 0, and M = Nm is the total mass 



In the present exploration, the force exponent spans the range 
—5/2 ^ a ^ 5/2, and the initial conditions are characterized 
by values of the virial ratio < r/o — 2Ki ti /\Wi a \ < 0.5. The 
N particles are distributed in space with a standard rejection 
technique. In the first family, used to study the evolution 



of cuspy initial conditions, we adopt the Hernquist ( 1990 1 
density profile 



p(r) 



Ma 



2irr(r + a) 



(16) 
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where M is the total mass, a — r*(v2— 1) is the scale radius, 
and r* is the half mass radius. In the second family of initial 
conditions, characterized by a flat core, we use the |Plummer| 
(19111 density profile 



p(r) = 



3Ma 2 



4tt (r 2 +a 2 ) 5/2 ' 



(17) 



in this case a = r»%/2 2/ ' 3 — 1. We then extract the initial 
velocity of the particles from a position-independent Gaussian 
distribution, with the velocity dispersion tuned as to obtain 
the desired initial virial ratio 770. The results are independent 
on the values of M and a, which do not appear in the 
dimensionless equations of motion (equation 12), so for each 
of the two families we explore the two dimensional parameter 
space defined by the pair (0,770)- We note that similar initial 
conditions have been used recently for the numerical study 



of violent relaxation in Newtonian gravity (e.g. |Visbal et al.| 
|2012| and |Sylos Labini|2013[ ). All the simulations presented in 
this paper were performed on a cluster of LINUX HP@Z700 
workstations, and each run (on a single processor) lasts for 
« 4 days when extended up to 50 t*. 



3.3 Analysis of the numerical outputs 

The numerical outputs are used not only to study the depen- 
dence of the virialization process on the force law but also 
to investigate the structural and dynamical properties of the 
final states. We assume that the system has reached its final 
state when the amplitude of the oscillations of 77 becomes 
smaller than 10 -4 (which typically occurs at t « 30/;*). Fol- 



lowing jMpotTeFaL] (|2002J) and |Meza fc Zamorano] (1997 1 
we compute the second order tensoij^ 



r ( fe ) r 



(k) 



(18) 



for the particles inside the sphere of radius rss containing 
the 85% of the total mass of the system, where r* are 
the Cartesian components of the position vector in the 
reference frame with origin in the centre of mass. The 
matrix 7y is diagonalized iteratively, requiring that the 
percentage difference of the largest eigenvalue between 
two iterations to be smaller than 10 -3 . This procedure 
requires on average 10 iterations, and we call 7i ^ I2 ^ I3 
the three eigenvalues. We finally apply a rotation to the 
system in order to have the three eigenvectors oriented 
along the coordinate axes. For of a heterogeneous ellipsoid 
of semiaxes a, b and c, we would obtain I\ = Aa 2 , I2 = Ab 2 
and 73 = Ac 2 , where A is a constant depending on the 
density profile. Consistently, for the end-products we define 
b/a = \fhJT\ and c/a = yflsJTi, so that the ellipticities in 
the principal planes are ei = l — y/h/Ii and €2 = 

It is known (see e.g. |van Alba da 1982, |Londrillo et al. 



1991 Trenti et al. 2005 N06a) that the end products of 



Newtonian dissipationless collapses starting from cold initial 
conditions (and also of MOND dissipationless collapses, see 



3 Notice that lij is not the inertia tensor, which is given instead 
by Trih^Sij - 



N07a) have surface density profiles well described by the 
Sersic law 



E(J?) = E e e 



(19) 



where b ~ 2m - 1/3 + 4/405m (Ciotti fc Bertin||l999| >, and 
E e is the projected mass density at effective radius R e , the 
radius of the circle containing half of the projected mass. In 
practice, in our analysis we circularize the projected density 
in the 3 principal planes of the virialized systems, and by 
particle count we determine the corresponding pair (R e , E e ), 
so that from eq. ( 19 I we obtain the best-fit m. In this way 
for each simulation we determine 3 sets of (E e , R e , m) and 
we then chose randomly one of them, being the others in 
general qualitatively similar. 

In the same spirit as previous works (DC09, DCC11) 
we focus on different indicators of relaxation, such as the 
evolution of the virial ratio and of the phase-space sections 
(r, v r ). We also study the final differential energy distribution 
n(E), defined by the relation 



i(E)dE = N 



(20) 



(e.g. Binney & Tremainc 2008). Finally, we construct the 



so-called pseudo-phase-space density of the final states and 
we check if they obey the GDSAI, as described in Sect. 4.3. 



3.4 Testing the code 

As a first set of numerical experiments, we determined the 
optimum choice of the softening length and of the time 
step to be used in the simulations. Following the procedure 
described in Section 3.1, we found that, independently of 



10" 



guarantees not only numerical accuracy of 



the results (with energy conservation better than 3% at 
virialization in the worst cases, and usually better than 
the 0.5%), but also acceptable computational times. In 
addition, comparing the evolution of collapses starting from 
identical initial conditions with adaptive or fixed timesteps, 
we found that a fixed timestep At ~ 10 -2 t* guarantees a 
good balance between computational time and energy and 
total angular momentum conservation independently of 
the value of a and of the initial profile, so we adopt this 
criterion for all the simulations. 

We tested our direct code in the Newtonian case against 
several well established results of numerical simulations of 
cold collapses obtained with the tree-code FVFPS ( |Lon-| 
drillo et al.|2003[ ), as well as in the Newtonian limit with the 



parti cle-mesh MOND code N-MODY ( |Londrillo fc Nipoti 
|2009[ ). In particular, we performed collapses for different ini- 
tial density profiles and values of the virial ratio. We fit the 
final surface density profile with the Sersic law over the radial 
interval 0.37? e — 107? e . As shown in Fig. [I] (left panel) the 
resulting values of m range from ~ 2 for the hottest initial 
condition (r/o = 0.5) to ~ 4 for the cold collapse (r/o = 0). 
Over the radial range here considered, the final profiles are 
indistinguishable from those obtained by N06a for compara- 
ble values of the initial virial ratio (0 ^ 770 ^ 0.2). As can 
be seen from Fig. [2] the final states, for both Hernquist and 
Plummer initial profiles, are roughly spherical for 770 > 0.1 
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Figure 3. Time evolution of the virial ratio n for cold (r/o = 0) Hernquist initial conditions with different values of a. Upper panel: 
gravity-like forces (a > 0). Lower panel: harmonic oscillator-like forces (a ^ 0). For clarity, the unrelaxing case of the harmonic force 
(a = —1, solid line) has been plotted only up to the third peak of rj. 



and prolate (c/a ~ b/a — 0.5) for rjo ^ 0.1, consistent with 
the results of N06a, which indicates that colder systems are 
more prone to undergo instabilities that perturb significantly 
their initial shape. Finally, the tests confirm that after virial- 
ization n(E) is well described over a broad range of energies 
by an exponential function. Remarkably, slightly bimodal 
final differential energy distributions characterize the systems 
starting from hotter initial conditions (Fig. [I] right panel) as 
seen in analogous plots of previous papers (figure 2 in N06b 
and figure 8 in Londrillo et al.|1991 l. 



4 RESULTS 

4.1 The relaxation process 

One of the motivations of this study is to elucidate the reason 
of the long relaxation time (in units of their dynamical time) 
of dMOND systems when compared to the same quantity 
for Newtonian systems. A simple measure of the relaxation 
effectiveness can be obtained considering the number and 
the decay rate of the major oscillations of the virial ratio 
rj. In Fig. [3] we show the evolution of 77 as a function 
of the dimensionless time r for cold (rjo = 0) Hernquist 
initial conditions with different values of a. The top panel 
illustrates the evolution for the family of gravity-like forces. 
The similarity with the results obtained with the shell 
models (Fig.l in DCCU) is remarkable: a decrease of a 
leads to a higher value of rj at the first peak, and to a 
longer series of virial oscillations of decreasing amplitude. 



Curiously, for a = 1 the amplitude of the virial oscillations 
increases again at large times, similarly to dMOND collapses 
(N07a, |Ciotti et al.||2007[ ), confirming that dMOND behaves 
qualitatively as the 1/r force when considering a system 
not deviating too much from spherical shape, reinforcing 
the previous result of the long relaxation times of N— body 
systems governed by the non-linear field equation of MOND. 
We interpret the large peak values of rj for small a as due 
to the fact that the force is stronger on large scales, and 
that the systems with low a collapse more as a whole, 
consistently with the force being more similar to the 
harmonic oscillator case. It is important to recall that in the 
Newtonian case a spherical homogeneous shell does not exert 
any force inside, while inside a shell the force is directed 
outwards for a > 2, and the opposite happens for a < 2 (e.g. 
DCCU). Therefore, when a > 2 the external regions act 
against the collapse, while for a < 2 also the external re- 
gions of the system contribute more and more to the collapse. 

The results for collapses driven by harmonic-like forces 
are shown in the bottom panel of Fig. [3] Note how the virial 
ratio oscillates with peak values of rj significantly larger, and 
more regular oscillations than in the gravity-like cases. As 
expected, in the a = — 1 case no relaxation takes place, since 
the whole system behaves as a single harmonic oscillator 
(e.g. |Lynden-Bell fc Lynden-Bell|[ll82l see also eq. (12}). 
In particular, in a system of harmonic oscillators starting 
at rest all the particles cross the centre simultaneously, so 
that \W\ — > while K — > —E, causing rj to diverge. For 
the reasons described above, in the super-harmonic case 
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Figure 4. The evolution of phase-space sections (r,v r ) for Hernquist initial conditions with t]q = 0, and —2 ^ a ^ 2. As expected, no 
mixing is acting on the harmonic-oscillator case (a = — 1, where the straight line rotates clockwise and each point of it describes similar 
ellipses), while mixing appears again in the superharmonic case (a = —2). The characteristic time and length scales, t» and r*, are related 
as given by eq. (10). 
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Figure 5. Left: projected density profiles of the end-products of Hernquist initial conditions with rjo = and different values of a (dots). 
Dashed lines are their Sersic best fit. Right: the associated pseudo phase-space density. Residuals with respect to the best fit are also 
shown. 



(a < —1), the first peak of 77 is reached at earlier times 
for decreasing a. However, the peak values decrease, due 
to phase mixing which forbids all the particles to cross 
the centre simultaneously. Cases with rjo slightly larger 
than 0, and for the Plummer initial conditions, are not 
shown, being qualitatively the same. In general, large val- 
ues of r/o correspond to small amplitudes of the first peak of r\. 



The long lasting virial oscillations for forces close to the 
harmonic oscillator are associated with a poorer mixing in 
phase space. Such behavior is evident from the evolution in 
the phase-space section (r, v r ) defined as radial position and 
radial velocity. In Fig. [4] we show snapshots of the phase- 
space at 1, 10 and 50 t* for a —-2, -1, 1, 2. Consistently 
with the findings of DCC11 (where the narrower interval 
of 1 ^ a sC 2 was studied), larger values of a show a more 
efficient phase mixing with respect to systems with a ^ 0. 
Again, the similarity with the plots in N07a (their figure 4) 
and Ciotti et al. (2007, their figures 2 and 3) is remarkable. 
One may speculate that the coherent structures in phase 
space that persist at large times (in units of £*) are akin 
to the so-called phase-space holes reported by some authors 



(Mineau et al. 1990 Joyce & Worrakitpoonpon 2011 and 
|Teles et al.||2011[ ) in the context of the one dimensional 
infinite sheet model, where the mixing is quite poor. It must 
be pointed out that the poorer mixing in Newtonian gravity 
in lower dimensions is essentially due to the smaller number 
of degrees of freedom that are involved in the relaxation 
rather than a different exponent in the force law (see e.g. 
|Kandrup|19 89 ). As already remarked, the efficiency of phase 
mixing is non monotonic with a, with no mixing for a = — 1. 



4.2 Structural properties of the end products 

The triaxiality of the final states of the collapses is shown 
in Fig. [2] where we plot the values of the axial ratios b/a 
and c/a of the end products at 50 i„, for representative 
values of a and for increasing values of 770. In general, 
the Newtonian behavior is confirmed, in the sense that 
at fixed a triaxiality is more pronounced for small values 
of 170, for both Plummer and Hernquist initial conditions. 
Again, the only exception is the a = — 1 force, when the 
systems retain their spherical shapes, consistently with 
their orbital structure. For given r/o, the triaxiality as a 
function of a shows characteristic non-monotonic trend 
especially visible in the bottom panel of Fig.[2]for rjo = and 
770 = 0.1. For decreasing a the flattening increases, reaches 
a maximum, and then decreases again. The maximum 
values of the triaxiality (the minimum values of b/a and 
c/a) are obtained for 1 ^ a 2 when ^ rjo ^ 0.1, with 
a quite clear correlation between a and rjo- Remarkably, 
no models are found to be flatter than an elliptical galaxy 
E7 (i.e. c/a ^ 0.3), thus leading to conjecture that this 
limit may hold more generally than just in Newtonian gravity. 

In analogy with the case of Newtonian collapses, we 
fitted the final projected density profile with the Sersic 
law as described in Sect. 3.3. In general, we find that the 
Sersic law provides a good description of most of the final 
states for both Plummer and Hernquist initial conditions. 
More quantitatively, in the left panel of Fig. [5] we show 
the end state of the collapse of a perfectly cold Hernquist 
initial condition. The main result is a quite well defined 
dependence of the Sersic index m on a, with large values of 
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Figure 6. Projected density profiles of the end products at 50 t* (points) and their Sersic best fits (lines) of Plummer initial conditions 
for a = — 2, —1, f , 2.5 and different values of rjo. In the a = — 1 case (harmonic force), the projected density profile is scale invariant, as 
expected. The symbols are the same as in Fig. 1. 



m associated with large values of a. In practice, gravity-like 
forces produce more peaked density profiles than harmonic 
like forces. The a = — 1 case is not shown, as the profile 
collapses and expands self-similarly. For other values of 
a, percentual deviations of the data from the fits are well 
within the 20%. The model with the largest deviations is 
the superharmonic one with a = —2 and m ~ 2.4. The 
dependence of the final states on 770 is shown in Fig. [6] In 
general, hotter initial conditions lead to smaller values of m, 
independently of a (for the Newtonian case see Fig. [l] left 
panel). Moreover, the largest deviations from the best fit are 
again produced by the superharmonic a — —2 force, while 
the a — 1 case is remarkably similar to the dMOND results 
of N07a. 

For completeness, in Fig.[7|we also show the three dimen- 
sional (angle-averaged) density profiles of the end-products 
of cold (770 = 0) Plummer initial conditions. As apparent, 
and in agreement with the projected density profiles, higher 
valiues of alpha corresponds to more peaked final density 
profiles, while for a = — 1 the normalized profile does not 
change. 



4.3 The differential energy distribution, the 

pseudo phase-space density and the GDSAI of 
the end products 

In addition to their structural properties, the virialized 
final states of collapses are usually also studied from the 
point of view of the phase-space properties. Here, following 
a well established approach, we focus on their differential 
energy distribution, on the radial trend of the so-called 
pseudo phase-space density, and finally on the density-slope 
inequality. 

In Fig.[| we show the final differential energy distribution 
n(E) for —2 ^ a ^ 2.5, and for different values of 770- Each 
distribution is normalized to the total number of particles, 



T 




Log(r/r h ) 

Figure 7. Normalized angle averaged three dimensional density 
profiles of the end products of Plummer initial conditions (black 
solid line) with rjo = 0, and different values of a. Radii are 
normalized to the volumetric half-mass radius of the final state. 

and the energy range to -B max (for a ^ 1) and to |-E m in| (for 
a > I) of the final states. The initial conditions with 770 = 
are represented by the heavy solid lines. The first important 
and general feature is that n(E) is peaked at high energies for 
a > I (see also Fig. 1, right panel, for the Newtonian case). 
In practice, the virialized final states of systems with gravity- 




Figure 8. Normalized differential energy distribution n t (E) = n(E)/N of the end products of Plummer initial conditions for different 
values of tjq and a = —2, —1, 1, 2.5. The reference energy is i? ma x or |£ m i n | if a ^ 1 or a > 1 respectively. The Newtonian case (a = 2) is 
shown in Fig.l. The heavy solid lines represent the n t (E) for the initial condition with r]o = 0. Hernquist initial conditions lead to very 
similar energy distributions. 



like forces allowing for escape are mainly supported by loosely 
bound particles (e.g., Binney & Tremainc 2008 Binney|1982 



Ciotti|1991| ) corresponding to particles in the outer regions 
Note also how n(E) evolves significantly due to relaxation, 
with major changes at high energies. For harmonic-like forces 
the situation is different, and very little evolution is found. 
Both the initial conditions and the final states have a n(E) 
distribution peaked at low energies. Of course, consistently 
with the extraordinary nature of the harmonic force, the 
n(E) for the a = — 1 case is not evolving (barring numerical 
fluctuations). In general, different values of r/o in the range 
explored do not affect significantly the shape of the n(E) with 
the exception of a = 1 forces. The systems with a = 1 show 
an intermediate behavior, with a trend similar to dMOND 
collapses, (See N07a, Fig. 5 therein). For this latter case it 
is apparent how decreasing values of r/o tend to populate 
the external regions of the final systems. As discussed in 
the Introduction, a specific feature of the n(E) obtained 
in Newtonian collapses is the exponential shape over some 
energy range. Here, due to the energy sign associated with 
the value of a (see discussion in Sect. 2) we consider the 
function 



n(E) = A x 



-b\e\ 

-BE 



a > 1; 
a € 1; 



(21) 



where and A are an inverse (positive) temperature and 
a normalization factor respectively. It is apparent that 
for a < 1 the shape of n(E) cannot be described by a 
single-temperature exponential distribution, independently 
of the hotness of the initial conditions. Instead, for 
gravity-like forces with a ^ 1 a larger energy range exists 
over which n(E) can be qualitatively described with an 



exponential function as in eq. (21 1. The main difference in 



the gravity-like forces is between a = 1 case and the other 
case with a > 1 (see also Fig. 1): while in the forces allowing 
for escape (a > 1) the exponential region is peaked towards 
high energies, in the a — 1 case the peak is at low energies, 
corresponding to the central regions. Interestingly, for a ^ 1 



the trend between the Sersic index m and the inverse 
temperature 6 of the best fit n(E), is qualitatively similar to 
what found by |Ciotti| ( |1991[ ) in the analysis of the Newtonian 
Sersic models. We finally note how the current N— body 
simulations produced final n(E) much better described 
by an exponential distribution than in the shell model 
(DCC11), a natural consequence of a better energy exchange 
among the components of the system. For the final states of 
the systems with a > 1 we also considered the fraction of 



escapers (i.e. particles having positive energy, see also Joyce 



et al.||2009| and Sylos Labini 20131. As expected, at fixed 



a and for given initial density profile, with low values of 
rjo (i.e., cold initial conditions) there is a larger number of 
escapers. Also, for fixed density profile and r/o, the fraction 
of escapers is found to be weakly dependent on the value of 
a, with values ~ 3% for most cases and with the maximum 



value of : 



I for perfectly cold Plummer model with a = 2.5. 



Another property of interest, recently focus of several 
investigations, is the so-called pseudo phase space density 
e.g. see Taylor fc Navarro|2001 Ascasibar fc Binney|200^ 



Hansen et al.|2010[ |Ludlow et al.|2010[|Barber et al.|2012 



Sparre fc Hansen|2012 I defined as 

Pjr) 
cr 3 (r) ' 



(22) 



where p(r) and er(r) are the angle-averaged density and 
velocity dispersion at radius r. For Newtonian collapses, the 
numerical simulations have unequivocally shown that Q is 
described quite well by a power-law 



X k 1.87. 



(23) 



neither the origin of this relationship nor the dependence 
of this property on the initial conditions are yet fully 
understood despite the numerous efforts. On one side, the 
results of numerical simulations seem to point out to a 
remarkable robustness of eq.(|23[): even though the power-law 
trend of Q was initially considered a peculiarity of the NFW 
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profiles ( Navarro et al.|1997 l, it is known that other profiles 
share this property (e.g. the family of self-consistent /oo 
models, see Bertin &: Stiavelli|1984 Zocchi|2010 1. However, 
there are self-consistent equilibrium systems where Q is 
not a power-law (e.g. the Plummer sphere). Therefore it is 
natural to ask wether the power-law is a specific feature of 
violent relaxation in Newtonian gravity or its origin should 
be searched more in the physics of dissipationless collapse, 
independently of the force law involved. Here we are in 
the ideal position to address this question, and in fact the 
obtained results are quite significant. As can be seen in 
Fig. [5] (right panel) a power-law trend for Q is reproduced 
surprisingly well also in the case of non-Newtonian forces. 
For all the considered case (with the exception of a = —1), 
and independently of the initial density profile, at fixed 
rjo the exponent \ increases for decreasing a and the 
function Q steepens. We also found that at fixed a, Q 
steepens for decreasing r]Q. These findings lead to conclude 
that a power-law radial dependence of Q is more a con- 
sequence of violent relaxation than of the Newton gravity law. 

Finally, we check if the (angle-averaged) Global Den- 
sity Slope Anisotropy Inequality (GDSAI) is obeyed by the 
end products for different values of a. Ciotti & Morganti 
(2010ab), prompted by the important asymptotic result of 
An & Evans ( |2006[ ), proved that a very large class of Newto- 
nian stellar systems with positive phase-space distribution 
function, necessarily obey the inequality 



7(r) > 20(r), Vr 



where 



7 = 



dlnp 



d In r 

is the logarithmic density slope and 



1 



A. 

2cr2 



(24) 



(25) 



(26) 



is the usual anisotropy parameter ( Binney & Tremaine^ 2008 ) . 
In the formula above oy and at are the radial and tangential 
component of the velocity dispersion tensor respectively. In 



particular, Ciotti & Morganti (2010b I speculated about a 



possible universality the GDSAI, even though their analytical 
methods where unable to treat the cases with /3(0) > 1/2. 
Significant progress and clarification has been made in the 
subject ( |van Hese et al.||2011 An et aL]|2012 1 and now the 
case of systems with separable augmented density is well 
understood: the GDSAI is obeyed by all separable systems 
with ,3(0) ^ 1/2, while counterexamples exist for systems 
with /3(0) > 1/2. Much less is known about systems with 
non-separable augmented density, but numerical simulations 
in Newtonian gravity seem to suggest that also in general 
systems the GDSAI is usually satisfied. Here we analyzed 
the results of the simulations for different values of a and rjo . 
The trends of 7 and /3, obtained using spherical averages and 
excluding the innermost regions (where discreteness effects 
dominate), revealed that the final states obey the GDSAI 
(with the obvious exception of the harmonic oscillator force). 
This is shown in Fig. [9] where, even in presence of numerical 
noise, it is apparent that overall 7 ^ 2/3, reinforcing the 
idea that the physical process leading to the establishment 
of the GDSAI may be independent of the specific force law 
considered. Finally, from the bottom panel of Fig. [9] it is also 



interesting to note how the final systems are significantly 
radially anisotropic in their outer regions, and become more 
and more isotropic near the centre, as also commonly found 
in Newtonian simulations. The a = — 2 case stands out as 
the less anisotropic and it is curious to recall that these 
systems are also those for which the Sersic law provides the 
less satisfactory description (Fig. |6j. 



5 DISCUSSION AND CONCLUSIONS 

As discussed in the Introduction, several theoretical 
arguments point toward the importance of elucidating 
the process of dissipationless collapse and virialization 
of iV— body systems with additive interparticle forces 
proportional to r~ a , a generalization of the Newtonian force. 
For this task we built a direct N— body code: preliminary 
results obtained with a shell model in spherical symmetry 
( Pi Cintio||2009| |Di Cintio fc Ciotti|r2011[ ) appear to be 
confirmed by the present simulations. The main results can 
be summarized as follows. 



The relaxation process, independently of the initial 
density profile (Hernquist or Plummer), is characterized by 
a first phase of strong oscillations of the virial ratio, followed 
by a gentler phase of relaxation. For decreasing a, the peak 
value of the virial ratio increases reaching a value formally 
infinite in the case of a perfectly cold collapse with a = — 1 
(i.e. a system of harmonic oscillators), and then decreases 
again. This non-monotonic behavior is a consequence of 
the different degrees of phase mixing as a function of a. 
Qualitatively, this effect can be understood by considering 
the force field inside a shell of matter for different values of 
a. As expected, systems with a = —1 do not relax due to 
their extraordinary orbital structure in which each particle 



behaves as an isolated harmonic oscillator (e.g. Lynden-Bell 
fc Lynden-Bell|1982 l. 



With the obvious exception of the a = — 1 force, when 
the systems retain their initial spherical shape, the final 
states are triaxial. As a rule triaxiality increases for colder 
initial conditions, similarly to what happens in Newtonian 
collapses. However, for fixed initial virial ratio 770, triaxiality 
is not a monotonic function of a: for decreasing a the 
flattening increases, reaches a maximum and then decreases 
again. The value of a for which the triaxiality is maximum 
depends on r/o, but it is almost independent on the initial 
density profile. Remarkably, no models are found to be 
flatter than an elliptical galaxy of type E7, independently of 
a and rjo- 

In general the Sersic law provides a good description of 
most of the final states, with large Sersic index m associated 
with large values of a (i.e. for gravity- like forces), and 
with cold initial conditions. Hotter initial conditions and 
harmonic-like forces produce density profiles characterized by 
smaller m. Moreover the quality of the Sersic fit deteriorates 
for low values of a. 

The differential energy distribution n(E) of the final 
states shows two distinct behaviors, separated by the case 
a = 1. In particular, n(E) is well described over a large 
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Figure 9. Radial profiles of the angle-averaged velocity dispersion (top), GDSAI indicator (middle), and anisotropy parameter (bottom), 
for the final virialized states of cold (rjo = 0) Hernquist initial conditions, is the volumetric half-mass radius of the final states. Similar 
trends arc found also for larger values of r)o and for Plummer initial conditions. Note the peculiar off-center maximum of a in the 
supcrharmonic force case. 



range of energies by an exponential function peaked at high 
energies for gravity-like forces with a > 1. For a = 1, the 
final n(E) is also exponential, but now the distribution is 
peaked at low energies (as in MOND simulations, N07a). 
Finally, when a < 1, very little evolution is found, and n(E) 
remains peaked at low energies. Remarkably, for a ^ 1 the 
trend of the inverse temperature 8 with the Sersic index m 
is similar to what found for the Newtonian Sersic models, 
with 9 increasing for decreasing m (Ciotti 1991). 



We found that the pseudo phase-space density 
Q — p/a :i of the (angle-averaged) final states is described 
very well by a power law of radius r _x , over a large 
radial range. In general, Q steepens for decreasing a, 
while at fixed a it flattens for increasing values of tjo ■ 
In addition we also found that the GDSAI holds (well 
within numerical uncertainties) for all the virialized 
end-states, and the amount of radial anisotropy tends to 
be higher for a > (i.e. for gravity-like force) than for a < 0. 



(in units of the internal dynamical time) is longer for a — 1 
than for a = 2, and this is due to the force becoming more 
similar to the harmonic oscillator case (a = — 1) when the 
system oscillates forever. Future explorations, exploiting the 
similarity between 1/r and dMOND forces, will be focused 
on the study of radial orbit instability for 1 /r a forces, in the 



line of the MOND study of INipoti et al. (2011 1 
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Overall the main conclusion of the present study are that 
several structural and dynamical features of the virialized 
states of cold and dissipationless collapses are not restricted 
to the special nature of the Newton law, but they appear 
to be more a property of the long range forces. Among the 
radial r~ a forces, however, the gravity-like forces (a > 0) 
are those with the results more similar to the 1/r 2 force. In 
addition, we found that the systems with interparticle force 
proportional to 1/r behave in many respects as dMOND 
systems. In particular, we confirmed that the relaxation time 



References 

An J.H. & Evans N.W., 2006, ApJ, 701, 1500 
An J.H., van Hese E., Baes, M., 2012, MNRAS, 422, 652 
Andredakis Y. C, Peletier R. F., Balcells M., 1995, MNRAS, 
275, 874 

Ascasibar Y.& Binney J., 2005, MNRAS, 356, 872 
Barber J. A., Zhao H., Wu X., Hansen S. H., 2012, MNRAS, 
424, 1737 

Bertin G, Ciotti L., Del Principe M., 2002, A&A, 386, 149 



Relaxation of N -body systems with additive r a interparticle forces 13 



Bertin G. & Stiavelli M., 1984, A&A, 137, 26 
Bertin G. & Trenti M., 2003, ApJ, 584, 729 
Bekenstein J. & Milgrom M., 1984, ApJ, 286, 7 
Binney J., 1982, MNRAS, 200, 951 

Binney J. & Tremaine S., 2008 Galactic Dynamics, 2nd Ed. 

(Princeton University Press) 
Bouchet F., Gupta S.. Mukamel D., 2010, Physica A, 389, 

4389 

Brandao CS.S. & de Araujo J.C.N. , 2012, ApJ, 750, 29 
Caon N., Capaccioli M., D'Onofrio M., 1993, MNRAS, 265, 
1013 

Chavanis P.H., 2008, in Dynamics and Thermodynamics of 
systems with long-range interactions: Theory and Experi- 
ments, AIP Conf. Proc, 970, 39 

Ciotti L., 1991, A&A, 249, 99 

Ciotti L., 2000, Lecture Notes on Stellar Dynamics. Scuola 

Normale Superiore, Pisa 
Ciotti L., 2009, NCimR, 32, 1 
Ciotti L. & Bertin, G., 1999, A&A, 352, 447 
Ciotti L., Londrillo, P. & Nipoti, C, 2006, ApJ, 640, 741. 
Ciotti L. & Morganti L., 2010a, MNRAS, 401, 1091. 
Ciotti L. & Morganti L., 2010b, MNRAS, 408, 1070. 
Ciotti L., Nipoti C, Londrillo P., 2007, Proc. Int. Workshop 

on Collective phenomena in macroscopic systems, World 

Scientific 177. 

Courteau S., de Jong R.S., Broeils A.H., 1996, ApJ, 437, 21 
Dehnen W., 2001, MNRAS, 324, 273 

Di Cintio P.F., 2009, Master Thesis, Bologna University 
(DC09) 

Di Cintio P.F. & Ciotti L., 2011, IJBC, 21, 2279 (DCC11) 
Gabrielli A., Joyce M., Marcos B., 2010, Ph.Rev.Lett., 105, 
210602 

Gerhard O.E. & Spergel D.N., 1992, ApJ, 397, 38 
Graham A. & Colless M., 1997, MNRAS, 287, 221 
Graham A., 1998, MNRAS, 293, 933 
Henon M., 1964, Ann. d'Astroph., 27, 83 
Hernquist L., 1990, ApJ, 356, 359 

Hansen S.H. & Moore, B. 2006, New Astronomy, 11, 333 
Hansen S.H., Juncher D., Sparre M., 2010, ApJ, 718, 68 
Joyce M., Marcos B., Sylos Labini F., 2009, MNRAS, 397, 
775 

Joyce M. & Worrakitpoonpon T., 2011, Ph.Rev.E., 84, 1139 

Kandrup H.E., 1989, Ph.Rev.A., 40, 7265 

Lynden-Bell D., 1967, MNRAS, 136, 101 

Lynden-Bell D. & Lynden-Bell R. M., 1982, Proceedings: 

Mathematical, Physical and Engineering Sciences, Vol. 455, 

No. 1982, p. 475, The Royal Society 
Londrillo P., Messina A., Stiavelli, M 1991, MNRAS, 250, 

54 

Londrillo P., Nipoti C, Ciotti L., 2003, MSAIS, 1, 18 
Londrillo P. & Nipoti C, 2009, MSAIS, 13, 89 
Ludlow A.D., Navarro J.F., Springel V., Vogelsberger M., 
Wang J., White S.D.M., Jenkins A., Frenk C.S., 2010, 
MNRAS, 406, 137 
Malekjani M., Rahvar S., Haghi H., 2009, ApJ, 694, 1220 
Malekjani M., Haghi H, Jassur D.M.Z., 2012, New Astron- 
omy, 17, 149 

Marcos B., Gabrielli A., Joyce M., 2012, CEJPh, 10, 676 
Meza A., & Zamorano N., 1997, ApJ, 490, 136 
Milgrom M., 2010, MNRAS, 403, 886 
Mineau P., Feix M.R., Rouet, J.L., 1990, A&A, 228, 344 
Moffat J.W. & Sokolov I.Yu., 1996, Phys.Lett.B, 378, 59 



Navarro J.F., Frenk C.S., White S.D.M., 1997, ApJ, 490, 
493 

Nipoti C, Londrillo P., Ciotti L., 2002, MNRAS, 332, 901 
Nipoti C, Londrillo P., Ciotti, L., 2006, MNRAS, 370, 681 
(N06a) 

Nipoti, C, Londrillo, P. & Ciotti, L. 2006, Science and 

Supercomputing at CINECA, 122 (N06b) 
Nipoti C, Londrillo P., Ciotti L., 2007a, ApJ, 660, 256 

(N07a) 

Nipoti C, Londrillo P., Ciotti L., 2007b, MNRAS, 381, 107 
Nipoti C, Ciotti L., Londrillo P., 2011, MNRAS, 414, 3298 
Plummer H. L., 1911, MNRAS, 71, 460 
Prugniel P. & Simien F., 1997, A&A, 321, 111 
Sanders R.H., 1998, MNRAS, 296, 1009 
Sanders R.H., 2008, MNRAS, 386, 1588 
Sparre M. & Hansen, S.H., 2012, JCAP, 10, 49 
Srinivasan K., Mahawar H., Sarin V., 2005, ICCS, 

LNCS3514, p. 107-1 14, V.S. Sunderam et al. (Eds.) 

Springer- Verlag Berlin 
Stein E.M., 1970, Singular integrals and differentiability 

properties of functions, (Princeton University Press) 
Sylos Labini F., 2013, MNRAS, 429, 679 
Takizawa M. & Inagaki S., 1997, |arXiv:astro-ph/9702002| t T 
Taylor J.E. & Navarro J.F., 2001, ApJ, 563, 483 
Teles T.N., Levin Y., Pakter, R., 2011, MNRAS, 417, 21 
Trenti M. & Bertin G., 2005, A&A, 429, 161 
Trenti M., Bertin G., van Albada T.S., 2005, A&A 433, 57 
Trujillo I., Graham, A., Caon M., 2001, MNRAS, 326, 869 
van Albada T.S., 1982, MNRAS, 201, 939 
van Hese E., Baes M., Dejonghe H., 2011, ApJ, 726, 80 
Visbal E., Loeb A., Hernquist L., 2012, arXiv:astro- 

ph/1206.5852 

Youngkins V.P. & Miller B.N., 2000, Ph.Rev.E, 62, 4583 
Zocchi A., 2010, Master Thesis, Milano University 



